Boundary-value problems in cosmological dynamics 
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The dynamics of cosmological gravitating system is governed by the Euler and the Poisson equa- 
tions. Tiny fluctuations near the big bang singularity are amplified by gravitational instability 
into the observed structure today. Given the current distribution of galaxies and assuming initial 
homogeneity, dynamical reconstruction methods have been developed to derive the cosmic density 
and velocity fields back in time. The reconstruction method described here is based a least action 
principle formulation of the dynamics of collisionless particle (representing galaxies). Two observa- 
tional data sets will be considered. The first is the distribution of galaxies which is assumed to be 
an honest tracer of the mass density field of the dark matter. The second set is measurements of 
the peculiar velocities (deviations from pure Hubble flow) of galaxies. Given the first data set, the 
reconstruction method recovers the associated velocity field which can then be compared with the 
second data set. This comparison constrains the nature of the dark matter and the relation between 
mass and light in the Universe. 

PACS numbers: 95.35.+d 



I. INTRODUCTION 

Cosmology is concerned with observing and modeling 
the universe on large scales: from our own Milky Way, 
other galaxies, galaxy clusters, super clusters up to the 
largest scales as probed by measurements of the cosmic 
microwave background radiation (CMB). These observa- 
tions span a huge range of scales and all strongly suggest 
that: 1) the dominant form of matter is dark (a factor 
of 6 in mass over the normal baryonic matter), 2) the 
clustering amplitude decreases with scale, and 3) struc- 
ture forms by gravitational amplification of ting initial 
fluctuations. These are some of the main component of 
the standard paradigm in cosmology. Violation of any 
of them or all of them is consistent with only a very 
limited set of observations, if any. Cosmology has had 
a great impact on other fields of physics and science in 
general. The shear existence of the gravitationally domi- 
nant dark matter has stimulated scientists' (and others') 
vivid imagination for a few decades now. Abundance and 
masses of non-standard particles have been constrained 
from the observed clustering pattern alone. In addition 
to gravity, hydrodynamical processes can greatly influ- 
ence the formation and evolution of galaxies, groups and 
clusters of galaxies. Hydrodynamical effects, however, 
play a minor role in shaping the observed distribution 
of galaxies on scales a few times larger than the size of 
galaxy clusters. Therefore, gravitational instability the- 
ory directly relates the present-day large scale structure 
to the initial density field and provides the frame-work 
within which the observations are analyzed and inter- 
preted. Gravitational instability is a non-linear process. 
Analytic solutions exist only for configurations with spe- 
cial symmetry, and approximate tools are limited to mod- 
erate density contrasts. So, numerical methods are nec- 
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essary for a full understanding of the observed large scale 
structure of the universe. There are two complementary 
numerical approaches. The first approach relies on N- 
body techniques designed to solve an initial value prob- 
lem in which the evolution of a self-gravitating system 
of massive particles is determined by numerical integra- 
tion of the Newtonian differential equations. Combined 
with semi-analytic models of galaxy formation, A-body 
simulations have become an essential tool for comparing 
the predictions of cosmological models with the observed 
properties of galaxies. Because the exact initial condi- 
tions are unknown, comparisons between simulations and 
observations are mainly concerned with general statis- 
tical properties. The second approach aims at finding 
the past orbits of mass tracers (galaxies) from their ob- 
served present-day distribution. The orbits must be such 
that the initial spatial distribution is homogeneous. This 
approach is very useful for direct comparisons between 
different types of observations of the large scale struc- 
ture. Most common are the velocity-velocity (hereafter 
v-v) comparisons between the observed peculiar veloci- 
ties of galaxies and and the velocity field inferred from 
the galaxy distribution in rcdshift surveys. This types 
of analysis yield the cosmological mass density param- 
eter f2 m . Any systematic mismatch between the fields 
serves as an indication to the nature of galaxy formation 
and/or the origin of galaxy intrinsic scaling relations used 
to measure the distances, provided that errors in the cal- 
ibration have been properly corrected for. This second 
approach also allows to perform back-in-time reconstruc- 
tions of the density field on scales ~ 5 h _1 Mpc [3J. 

Finding the orbits that satisfy initial homogeneity and 
match the present-day distribution of mass tracers is a 
boundary value problem. This problem naturally lends 
itself to an application of Hamilton's variational principle 
where the orbits of the objects are found by searching for 
stationary variations of the action subject to the bound- 
ary conditions. The use of the Principle of Least Ac- 
tion in a cosmological frame-work has been pioneered by 
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Peebles (1989) @ and has long been restricted to small 
systems such as the Local Group Q and the Local Su- 
percluster H . Early applications to large galaxy redshift 
surveys have been hampered by the computational cost 
of handling the relatively large number of objects. Sub- 
sequent numerical applications speeded up the method 
and allowed the reconstruction of the orbits of ~ 10 3 
particles (Shaya, Peebles & Tully 1995). However, it was 
only recently that the improvement of the minimization 
techniques and the use of efficient gravity solvers made 
it possible to deal with more than 10 4 objects 5], compa- 
rable to the number of objects contained in the largest 
all-sky galaxy catalogs. 



II. COSMOLOGICAL DYNAMICS 

For the background cosmology we work with a 
Friedmann-Robertson- Walker Universe. In this uniform 
background, the physical distance, r, between two points 
is r oc a(t) where a{t) is the scale factor. We consider a 
matter dominated universe with mean density p = Clp c 
with p c = 3H 2 /8irG. For a CI = 1, we get a critical 
density flat universe with a ~ t 2 / 3, . The Universe is ge- 
ometrically open for CI < 1 and close for Q > 1. Cur- 
rent observations indicate that the Universe contain a 
cosmological constant which makes it flat even though 
Cl ~ 0.3[9]. Apart from the dependence of a on t the 
presence of a cosmological constant has very little effect 
on our description here. In particular, the equations of 
motion of perturbations remain correct. We further de- 
fine, H(t) = a/a is the Hubble function, and denote the 
comoving coordinate of a patch of matter by x — r/a. 
The fluctuations are described by the density contrast 
S(x, t) — p(x,t)/p(t) — 1 and the comoving velocity by 
v = dx/dt. Also, let D{t) be the linear density grow- 
ing mode normalized to unity at the present epoch, and 
f(fl m ) = dlnD/dlna w fl° n 6 (e.g., Peebles 1980). The 
equations governing the evolution of fluctuations in a col- 
lisionless mass component in an expanding Universe are, 
The Euler equation, 



dv 
dt 



2Hv + v- Vv = -Vcpg, 



the continuity, 



OS 
dt 



V- (1 + S)v = 



and the Poisson equation, 

V 2 tp g = inGpS . 



(1) 



(2) 



(3) 



The term 2Hv in the Euler equation is due to the expan- 
sion of the cosmological background. The source term 
in the Poisson equation represents density fluctuations 
above the mean background density. 



A. Linear gravitational instability 

Neglecting the non-linear terms v ■ Vv and V • 6v, the 
equations of motion reduce to 



1 



f(Q)H 



V • v 



and 



S + 2HS = -H 2 C15 , 



(4) 



(5) 



where an over-dot indicates a time derivative. For a criti- 
cal density Universe (Ct = 1 and H = 2/3t), the equation 
([5]) gives Si oc t 2 / 3 and S2 oc t -1 , as the growing and 
decaying solutions, respectively. A few things to note. 
First, without the term 2HS the solutions would be expo- 
nential functions rather than power laws in time. Second, 
even in the linear regime, the decaying mode prevents a 
full recovery of the initial conditions, at t ss near the 
big bang cosmological singularity. Indeed, recovering this 
mode requires a precise knowledge of the present S and 
S (or t>), in order to prevent a blow-up as t — ■> 0. 

The relation Q has a simple interpretation. Since 
H ~ 1/t and t V 2 g ~ —S, it gives the intuitive rela- 
tion v ~ — V(f> g t between the acceleration, — V^> 9 and 
velocity. The relation has played a prominent role in 
the analysis of large scale structure. The density con- 
trast S(x) as estimated form the distribution of galaxies, 
could be used in this relation to obtain the associated 
peculiar velocity v(x). This velocity fields could be com- 
pared with the actual observed velocities of galaxies. A 
good agreement between the fields yields the cosmologi- 
cal density parameter, fl, and also a confirmation of the 
gravitational instability mechanism for structure forma- 
tion. But, perhaps more interestingly, any mismatch be- 
tween the fields could be an indication of strange mode of 
galaxy/structure formation the result of which is a galaxy 
distribution different from that of the dark matter. 



B. Non-linear cosmological dynamics 

Linear theory is valid only when the fluctuations are 
small. In practice this is achieved by smoothing the ob- 
served galaxy distribution on small scales ( <10Mpc). 
We describe here some non-linear methods which can be 
used for a variety of purposes, e.g. recovery of the ini- 
tial conditions, estimating v from the galaxy distribution 
and constraining the masses of galactic halos. Here we 
focus on the estimation of v. One can use numerical sim- 
ulations of non-linear gravity to calibrate semi-analytical 
non-linear generalizations to (TjQ). The approach is use- 
ful as it provides partial differential equations which can 
be solved for v for a given source term, S. Nevertheless, 
such generalizations are usually statistical in nature. In 
the following, we will describe a more rigorous and accu- 
rate approach. 
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We switch to a Lagrangian description for a system of 
N equal mass particles in an expanding universe. Each 
particle represents a patch of matter which, for practical 
purposes, could be a galaxy. The equations of motion are 
(% = !■>. N), 



dt 



2Hv, 



where g = — V0 g and is given by 



G 



\X Xj\ 



-Gpax , 



(6) 



(7) 



The equations can be derived from the action, 



S = 



dt£ 



G 
a 



E 



2n 



Xi X j 



pa^xf 



(8) 



under stationary first variations of the orbits that leave 
x fixed at the present epoch and satisfy the constraint 
t 1 / 3 ^ — > const as t — > [3, @|. The second condition on 
the velocities guarantees homogeneity near the big bang 
singularity t — > 0, preventing a blow up of the solutions. 
We expand the orbits in a time dependent base functions 
q n (t) in the form, 



Xi^t^j — X- 



i,0 



(9) 



where flf^o is the position of the particle i at the present 
epoch, and the vectors Ci >n are the expansion coefficients 
with respect to which the action is varied, i.e., they sat- 
isfy dS/dC^ n = 0. The base functions q n are chosen such 
that the boundary conditions are satisfied. 

Our strategy is to find orbits that are as close as possi- 
ble to the Hubble flow. Therefore, we search for the min- 
imum of the action and do not look for stationary points 
which might describe oscillatory behavior of the orbits. 
To find the coefficients C^ n that minimize the action, 
we use the Conjugate Gradient Method (CGM) which is 
fast and easy to implement. The gravitational force g 
and its potential are computed using the TREECODE 
gravity solver. The time integration in the expression for 
the action is done using the Gaussian quadrature method 
with 10 points at the time abscissa. The CGM requires 
an initial guess for Ci >n . We will use the term FAM, 
for Fast Action Method, to refer to the reconstruction 
method described here. In the standard FAM applica- 
tion we compute the initial guess using the linear theory 
relation between the velocity and mass distribution. The 
minimum of the action proved to be rather insensitive to 
the choice of initial guess for Ci. n , as we have checked by 
running FAM experiments with initial Ci. n both set to 
zero and to random numbers with appropriate variance. 



Besides the initial set of Ci. n , the other free parameters 
are the softening used by the gravity solver and the tol- 
erance parameter that sets the convergence of the CGM 
method. The success of the least action reconstruction 
method is illustrated in figure (P). 



III. DISCUSSION 



The rapid rotation of galactic disks revealed the ex- 
istence of dark matter halos which engulf the luminous 
component. The measured virial motions of galaxies in 
clusters of galaxies also require the a gravitationally dom- 
inant dark component. Away from bound systems of 
galaxies and galaxy clusters, field galaxies show coher- 
ent flow pattern which deviates from a pure Hubble ex- 
pansion. This coherent velocity field is a direct probe 
of the large scale dark matter distribution in as much 
as rotational speeds and virial motions are a measure of 
the dark matter in galaxies and clusters. Indeed, the 
cosmic gravitational field responsible for the motions of 
galaxies, mainly depends on the gravitationally dominant 
mass density field of the dark matter. The actual distri- 
bution of galaxies may well be quite different from the 
dark matter distribution. Recent analysis of the galaxy 
surveys, however, reveal a good match between the statis- 
tical properties of the galaxy distribution and the corre- 
sponding properties for the dark matter as inferred from 
numerical simulations of dark matter evolution in the 
universe. This is encouraging, but there may still be sig- 
nificant deviations between the distribution of the dark 
and luminous components, which are not reflected in sta- 
tistical comparisons. The only way to detect such devia- 
tions is via direct detailed comparisons between the mea- 
sured velocities of galaxies and velocities estimated from 
the galaxy distribution. These comparisons have been 
done in the linear regime. The overall agreement between 
the fields is impressive, but minor persisting mismatch is 
detected in some regions in the local volume. It is possi- 
ble that non-linear analysis a al the least action principle 
could mitigate some of the disagreement. This remains 
to be seen. The least action principle could also be used 
to recover the initial conditions [2] , allowing us to answer 
one of the fundamental question of whether or not initial 
fluctuations were gaussian [4[. 

The program is not without flaws. Many physical ef- 
fects need to be addressed in detail. Most pressing is 
incorporating the assembly (or merging) history of galax- 
ies. Galaxies reside in dark matter halos which form in a 
hierarchical manner from small to large. Thus our own 
Milky Way galaxy for example, is likely to have had a 
major merging activity some 8 Gyr ago. All reconstruc- 
tion methods assume that galaxies are point tracers of 
the mass density field and do not account for merging 
effects. 
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FIG. 1: Maps of 2D-projected peculiar velocities for points residing in a slice of thickness 6 h~ Mpc cut through a simulated 
catalog. The length of the vectors is drawn in units of 1 h _1 Mpc = 50 Kms -1 . The top row shows the least action predicted 
velocities (labeled FAMz). iV-body velocities are shown in the middle row. The velocity residuals, VivbodH — Vfamz, are displayed 
on the bottom. The maps shown in the panels to the left hand side refer to all the points in the slice while only the velocities 
of points with moderate density contrast are plotted in the central and right columns. 
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